Localized basis sets for unbound electrons in nanoelectronics 
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It is shown how unbound electron wave functions can be expanded in a suitably chosen localized 
basis sets for any desired range of energies. In particular, we focus on the use of gaussian basis 
sets, commonly used in first-principles codes. The possible usefulness of these basis sets in a first- 
principles description of field emission or scanning tunneling microscopy at large bias is illustrated 
by studying a simpler related phenomenon: The lifetime of an electron in a H atom subjected to a 
strong electric field. 



I. INTRODUCTION 

The theoretical description of the field emission (FE) 
mechanism has been traditionally based on the Wentzel- 
Kramers-Brillouin (WKB) approximation and simple 
band-structure models for the emitting tips 1 . Research 
on field-emitting graphitic compounds has attracted 
much attention ever since its first observation in car- 
bon nanotubes (CNT's) by de Heer et alA Because of 
their high aspect ratio as well as mechanical and chemi- 
cal stability, carbon nanotubes (CNT's) are regarded as 
potential materials for field emitters and this has been 
experimentally demonstrated in recent years 3 ^^. The 
electronic structure of CNT's and other graphitic mate- 
rials strongly depends on details at the atomic level and 
not only on the overall shape. Therefore, an atomistic de- 
scription of the emission process in these materials is per- 
tinent. In this regard, tight-binding^, static density func- 
tional theory (DFT) 7 - 8 - 9 , or time-dependent density func- 
tional theory (TDDFT) 10,11 calculations have been able 
to explore the theoretical possibilities of some graphitic 
compounds as field emitters. 

In the last decade a great deal of work has also been 
devoted to master electronic transport at the nanoscale. 
From a computational viewpoint the development of ab 
initio codes for quantum transport calculations has revo- 
lutionized the field, since, for the first time, reliable pre- 
dictions are possible for the resistance of nanoscale ob- 
jects when a current is driven throug h 12 ' 13 i 14 ' 15 i 16 . Many 
of these codes are based on the non-equilibrium Green's 
function formalism (NEGF) which, in its most com- 
mon implementation, requires the use of localized basis 
sets. This is an important drawback which hampers the 
straightforward applicability of these codes to describe 
FE processes where electron wave functions are unbound 
and extend into the vacuum. 

A precise description of extended wave functions is also 
important when the electron transport takes place be- 
tween electrodes or metallic tips that are significantly 
separated from each other. This is typically the case 
in scanning tunneling microscopy experiments where the 
metallic tip is several A away from the surface and the 
electron current becomes a tunneling current. When one 
goes from the contact regime to the tunneling regime the 



vacuum region becomes so wide that diffuse functions 
centered on the tip and the surface atoms cannot describe 
this region properly. Care must be taken complementing 
the basis set away from the tip and substrate 1 ^ for typ- 
ical tip-surface distances. This is particularly important 
when a large voltage is applied between tip and sample 
which can affect the shape of the tunneling barrier to the 
point that electrons are emitted into the vacuum before 
reaching the substrate^. 

In this report we examine the possibilities and lim- 
itations behind the use of localized atomic basis sets 
for the numerical description of free or unbound elec- 
trons. We are interested in the further implementation 
of these basis sets within our first-principles quantum 
transport package ALACANT (Alicante Ab initio Com- 
putation Applied to Nanotransport) 13 ' 14 i 19 ' 20 which in- 
terfaces GAUSSIAN0321 or CRYSTALS 2 packages. This 
is why we focus on the use of non-orthogonal gaussian- 
type functions. In Sec. II we present a detailed analy- 
sis of the completeness of gaussian basis sets in one di- 
mension within a desired range of energies. In Sec. Ill 
we generalize the analysis to three-dimensional gaussian 
wave functions. In Sec. IV, as a simple illustration, we 
compute the lifetime of an electron in the 2s state of an H 
atom in the presence of an electric field using an appro- 
priate basis set and compare with analytical calculations. 
Section V presents the conclusions. 



II. BASIC THEORETICAL CONSIDERATIONS 
IN ONE DIMENSION 

We consider first a one-dimensional free electron sys- 
tem described by the kinetic energy Hamiltonian T x = 

, .2 

2^2- Here and in the following we use atomic units 
(A.U.) where the electron mass is set to one. The eigen- 
states are thus propagating waves which in the coordinate 
representation are given by ipk{x) cx exp[z/cx] with a real 
wavevector k and an energy dispersion given by 



m = r 



(1) 



We want to describe free propagating waves with a ba- 
sis set of localized atomic orbitals. As mentioned in the 
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introduction, we focus on the case of gaussian-type func- 
tions which are employed by many first-principles codes. 
In one dimension (ID) a gaussian-type function is given 
by: 



4>(x) = -\/2a/7r exp [— ax 2 ]. 



(2) 



We now place an infinite number of gaussian orbitals 
on a regular one-dimensional grid with lattice spacing 
a, and define the ket vector \4>n) corresponding to the 
gaussian wavefunction localized at grid point x n — na e x 
by (x\<p n ) — <f>(x — na). The overlap integral between two 
gaussian orbitals on the grid separated by a distance na 
is thus given by 



s(na) := (4>\(f> n ) = cxp[— a(na) 2 /2] 



(3) 



For the kinetic energy integral between gaussian orbitals 
separated by a distance na we have 



t(na) := (<j)\f x \^ n ) 



1 



W 2 |<M 



-(1 - a(na) 2 )exp[ 



*(na) 2 /2}. (4) 



In this discrete basis set, the kinetic energy Hamiltonian 
T x can be diagonalized by transforming to Bloch states 2 ^: 



\<fik) 



E 



exp[ifcna] \<f> n ) =4> T x \ip k ) = e(k)\(p k ) (5) 



Multiplying the eigenvalue equation by (<p\ from the left 
and resolving for e(k), the following dispersion relation is 
obtained: 



e(k) 



■\T x \ip k ) Er 



t exp[ikna] t(na) 



T,n=-oo exp[ifcna] s(na) 



(6) 



Obviously, the finer the grid, i.e. the smaller a, the 
better the dispersion relation must approach that of free 
electrons. In fact, in the limit of a — > 0, the expression 
([6]) can be summed to yield the exact dispersion rela- 
tion E(k) independent of the gaussian exponent a (see 
Appendix) : 



[ dx t(x) exp[ikx] 

lim e(fc) = — — —— 

a^o J dx s(x) exp[ikx\ 



k 2 

k Y ^E(k). 



(7) 



Of course, this is somewhat trivial in reality since on an 
infinitely fine grid any orbital must expand the entire 
Hilbert space. A text-book choice in this limit is a basis 
of delta functions (Dirac basis), which is obtained from 
the gaussian functions by taking a — > oo: 



lim \l — exp[— a(x — na) A 

a — >oo V 7T 



S(x 



However, for all practical purposes a finite grid param- 
eter a > must be set. Our aim now is to find an optimal 
a for a given a. Note, that the "atomistic" description al- 
ways yields a finite maximal k- vector -the Brillouin zone 
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FIG. 1: Energy dispersions e(k) for a = 1A.U. and for dif- 
ferent values of a. The inset shows the maximal k-value fc max 
calculated from the inflection point of the energy dispersion 
e(k) as explained in the text. Fhe points have been calculated 
from individual dispersion curves while the line is a linear fit 
femax(a) = n/a — c • a to these points; c ~ 0.33A.U.. 



(BZ) boundary fcez = it /a— in contrast to truly free elec- 
trons where k is unbound. Fig. [1] shows the calculated 
dispersion e(fc) for fixed lattice parameter (a = 1 A.U.) 
and for different values of a. We see that in general 
the approximation works very well for smaller k- values, 
i.e. near the BZ center, but becomes worse towards the 
BZ boundary where the dispersion relation flattens out. 
Furthermore, the approximation is better, i.e. valid for 
a bigger range of fc's, the smaller a, i.e. the more diffuse 
the Gaussian function. However, it is already remark- 
ably good for relatively large values of a. For instance, 
for a = 1, e(k) is a good approximation for fc-values up to 
70% of the theoretical upper limit fcez ■ This can also be 
seen from the effective mass m*(k) = (d 2 e/dk 2 )~ 1 calcu- 
lated numerically from the dispersion relation e(fc), and 
shown in Fig. [2j The smaller a, the better does the ef- 
fective mass m* approximate the constant mass of truely 
free electrons (m e = 1). 

This suggests that e(k) — > E(k) for a — > for k < fcez, 
so that the optimal exponent a would be zero. Of course, 
setting a = is numerically not feasible. In fact, in 
most quantum chemistry codes the gaussian exponents 
a cannot be chosen arbitrarily small for computational 
reasons. On the other hand, as can be seen from Fig. 
[TJ the approximation is already quite good for reason- 
able values of a if we limit the range of wavevectors k 
to some A: max smaller than the theoretical upper limit 
/jbz- Since the slope of the real free-electron dispersion 
(8) increases linearly with k, dE/dk = k, we define k max as 
the fc-vector where the slope of the approximate energy 
dispersion e(k) starts to decrease, i.e. at the inflection 

point of e(k): c? 2 e/dfc 2 (fc max ) = 0. Clearly, at this point 
the effective mass m* = (d 2 e/dk 2 )~ 1 becomes infinite, 
so that our approximation is not valid any longer i^l As 
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FIG. 2: Effective mass m*(k) = {(fe/dk 2 )' 1 for a = 1A.U. 
and for different values of a. 



can be seen from the inset in Fig. [TJ the thus calculated 
&max can be fitted very well by a linear function of the 
gaussian exponent a that approaches the BZ boundary 
&BZ = tt/o for a — > 0. 

How does k max depend on the lattice parameter a? 
Suppose that for a given lattice parameter ao and ex- 
ponent ao we obtain a maximal range A:^ ax . Thus we 
have: 

k 2 

£ao,aoOo) ~ y fOT 8,11 fc ° - fc max- 

If we now alter the lattice spacing a and simulaneously 
scale a as a(a) = ao(ao/a) 2 and k as k(a) — ko(ao/a) 
it is straightforward to show that the energy dispersion 
is essentially unchanged apart from an overall factor of 
(a /a) 2 : 

«..«(*)=(-) e 0o>ao (fc )«(-J f = y 



^0 < Caxi 

Therefore 



Thus e a Q (fc) = fc 2 /2 for all fc with (a/ao)k 
i.e. for all k with /c < (ao/a)/;^ = fc m; 
fcmax scales exactly as the theoretical limit fcez = 7r/a 
like 1/a with the lattice parameter a if we simultaneosly 
scale the exponent a as 1/a 2 : k ma , x {a) — (ao/a)fcJ^ ax . 



III. THREE-DIMENSIONAL GAUSSIAN 
FUNCTIONS 



by 



In three dimensions a gaussian wavefunction is given 



$(r) 



(2a/7r) 3/2 exp[-ar 2 ] 



4>{x)4>(y)tf>{z). (9) 



Again, we place the three-dimensional (3D) gaussian 
functions on a regular ID grid along the x-axis with 
lattice spacing a, and define the ket vector |$ n ) corre- 
sponding to the gaussian function localized at grid point 
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FIG. 3: DOS projected onto a bulk site % (PDOS) of a finite 
(N=200I) chain of gaussian functions for different values of 
the gaussian exponent, a, together with the DOS for free 
electrons moving in one dimension. 



r„ = nae x by: 

<^$„) = $(r - r n ) = 4>{x - na)(f>(y)(t)(z) (10) 

Now there is an additional contribution to the kinetic 
energy integral from the two directions (y and z) per- 
pendicular to the direction (x) of the ID grid: 

t m (na) := ($\(f x + f y + f,)|$ n ) 

= (0|T|0 na )(0|0) 2 +2 x (0|0 na )<^)(0|f |$ 

= t(na) + 2 x s(na) x t(0) 

= t(na) + a s(na). (11) 

Then, using Bloch's theorem we obtain the corresponding 
dispersion relation for 3D gaussian functions which differs 
from the analogous expression for ID gaussian functions 
([Fj]) only by a constant energy shift: 



C3D(fc) = a + 



Er^L-oo t ( na ) exp[ifcna] 



y c 



s(na) exp[ikna] 



(12) 



The constant energy shift a is due to the lateral con- 
finement of the electrons to the region defined by the 
gaussian-shaped wavef unctions. Obviously, this confine- 
ment is eliminated by letting a — *■ 0, i.e. by increasing 
the diffuseness of the gaussian orbitals and thus letting 
the lateral wavefunction 4>(y) x <fi(z) become a free elec- 
tron wave with wavevector zero (k y = k z =0). However, 
as said before, in most quantum chemistry packages the 
gaussian exponents cannot become arbitrarily small for 
computational reasons. Thus when the artificial offset 
needs to be reduced beyond the computational limit for 
the exponents, one has to reduce the lateral confinement 
of the electrons by extending the vacuum basis set with 
gaussian wavefunctions along the other two dimensions. 

To analyze the density of states (DOS), we define the 
retarded Green's function in our non-orthogonal basis set 
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as (for details see e.g refs. |2Q||2 



G(E) = [{E + iS)S-T}-\ 



(13) 



where S is an infinitesimal quantity, S is the overlap ma- 
trix, and T is the 3D kinetic energy Hamiltonian rep- 
resented in the non-orthogonal gaussian basis set. The 
integral over the total DOS, p(E), must be equal to the 
number of basis functions, N, in our system. 



N 



dEp(E) 



1. 



-Im 



dETr[G(E)S] (14) 



In Fig. [3]we plot the DOS of a finite chain of N = 2001 
gaussian functions projected (PDOS) on a bulk site i, i.e. 
— ^lm[G(E) S]i^, for four values of the gaussian expo- 
nent a and a fixed grid parameter of a = 2 A.U.. The 
Green's function for an infinite system can be computed 
without difficulty 20 . We prefer, nevertheless, to present 
results for finite systems here since these are the ones 
employed in the next section. Obviously, no differences 
should be expected for large enough systems. As dis- 
cussed above, there is an unwanted contribution to the 
kinetic energy coming from the confinement in the two 
directions perpendicular to the direction of the atomic 
chain. This opens a gap in the PDOS close to zero en- 
ergy which disappears as the exponent, and thus the lat- 
eral confinement, decreases. For values of the exponent 
0.05 < a < 0.1 the PDOS already resembles that of free 
ID electrons for energy values from 0.0 to 1.25 A.U.. The 
small but still visible gap in the PDOS at zero energy due 
to the residual lateral confinement can be eliminated by 
letting a — > 0, but, as already mentioned, numerical lim- 
itations inherent to most computational packages do not 
recommend to do this. This already mentioned possibil- 
ity of increasing the cross section of the chain by adding 
functions laterally will be explored in the future. 



IV. IONIZATION OF THE 2s STATE OF THE 
HYDROGEN ATOM 

In order to test the use of gaussian basis sets in realistic 
FE calculations, we have studied the ionization probabil- 
ity of the first excited state, 2s, of the hydrogen atom 
in the presence of an electric field. For simplicity's sake, 
we have not included the 2p orbitals, although hybridiza- 
tion with them should be taken into account in a more 
realistic description of this problem. The calculations 
have been performed using the GAUSSIAN03 ab initio 
package as a complementary test for situations when the 
use of standard ab initio packages is convenient or neces- 
sary. The 2s orbital is coupled to two semi-infinite chains 
of gaussian functions representing the vacuum along the 
direction of the applied electric field as shown in Fig. 
H The straightforward use of GAUSSIAN03 forces us to 
work with finite vacuum chains as the one analyzed in the 
previous section. The results presented below have been 



FIG. 4: One dimensional model where a central 2s orbital of a 
hydrogen atom is connected to semi-infinite chains of gaussian 
functions on both sides representing the vacuum. The electric 
field is applied along the chain direction. 
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FIG. 5: Density of states projected on the 2s state of the 
hydrogen atom for different values of the applied electric field. 
In the inset the field dependence of the peak width is depicted 
and compared with the theoretical one (dotted line). 



converged in the length of the lateral chains (N — 1000) 
so that no finite size effects can be appreciated. 

The lifetime of the 2s state in the presence of an electric 
field can be extracted from the PDOS, 



P2s{E) 



-Im[G(E)S} 2 s,2s 



(15) 



The grid parameter and exponent for the two chains have 
been set to 1.0 A (~ 2 A.U.) and 0.07 A.U., respec- 
tively. The distance between the hydrogen atom and 
both vacuum chains has been fixed to 5.0 A (~ 10 A.U.) 
to optimize the hybridization between the vacuum states 
and the 2s state of the hydrogen atom. The 2s atomic 
orbital is modeled by the corresponding STO-6G basis 
function implemented in GAUSSIAN03 and the vacuum 
orbitals were implemented using ghost atoms. Care has 
been taken to remove the electron in the GAUSSIAN03 
input file, reducing thus the calculation to a simple non- 
interacting problem. 

The resonance width of the 2s orbital has been com- 
puted (see Fig. [5]) and compared to the theoretical one^S 
for some values of the electric field F (see inset in Fig[5]). 
A finite value of 8 (or fictitious width) is included to 
smear out the computed PDOS. This is why our results 
for the resonance width do not reach zero for zero applied 
field (see inset in FigJS]). As it would be expected the 
computed resonance width T2 S for the 2s state of the hy- 
drogen atom is well approximated by the quasi-classical 
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theory: 



r 2s (F) 



32F 2 



cxp 



12F 



(16) 



Considering the simplifications in our calculation which 
has been restricted to one dimension and that little care 
has been taken in completing the basis set around the 2s 
orbital, the results are qualitatively similar to the the- 
oretical ones. This illustrates that localized basis sets, 
when appropriately chosen, can be used to represent un- 
bound electron wavefunctions which are required to de- 
scribe FE phenomena or high-bias STM. 



V. CONCLUSIONS 

In this work we have shown the possibility of represent- 
ing unbound electron states using localized gaussian-type 
functions as a basis set. Although plane waves are a more 
natural basis for unbound electrons, localized basis sets 
have to be employed for the computation of field emission 
phenomena in the framework of many commonly used ab 
initio quantum transport packages such as ALACANT. 
In addition these basis sets provide a more natural way of 
quantifying local atomic properties when needed. As an 
illustration we have studied numerically the lifetime of an 
electron in the 2s orbital of the H atom in the presence 
of an electric field with the help of the GAUSSIAN03 
package. 
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APPENDIX A: DERIVATION OF EQ. © 

In the limit of small lattice spacing a the sums over the 
jrid points n in eq. ([6]) can be approximated by integrals 



over the quasi continuous variable x = na: 

oo If. 

} t(na) exp[ikna] « - / dxt(x)e lkx (Al) 



and 



°° If 

s(na) exp[ifcna] w - / dx s(x) e lkx . (A2) 

n— — oo ^ 

Thus the two sums are approximately given by the 
Fourier transform of the hopping integral t{x) and the 
overlap integral s(x), respectively. Since the overlap inte- 
gral is a Gaussian function, s{x) = exp[— ax 2 ], its Fourier 
transform s(k) is also simply a Gaussian: 

S(Jfe) = f dxs(x)e lkx = a~ 1/2 exp[~k 2 /2a}. (A3) 

V 2vr J 



The Fourier transform t(k) of the hopping integral 
5' 1 " >; 



t(x) = ~(1 — ax 2 ) exp[— ax 2 } is a bit more involved: 
1 



t{k) 



I2n 



dx t(x) e lkx 



dx x 2 s[x) e l 



2L S (k) -SLJ- 

2 W 2 v^F 
^-exp[-fc 2 /2a]-z 2 T ^ S (fc) 



-a 



-1/2 



1 k< 



ya- 1 / 2 exp[-fc 2 /2a]. 



exp[-fc 2 /2a] 



(A4) 



Thus, in total we obtain for the energy dispersion e(k) 
in the limit of small a: 



= J2 n t ( na ) expjikna] ^ i{k)_ = k? 
J2 n s(na) exp[ifcna] s(k) 2 

This proofs eq. ([7]) since in the limit a — > the approxi- 
mation of the sums by integrals becomes exact. 
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